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While the basic principles of conventional solar cells are well understood, little attention has 
gone toward maximizing the efficiency of photovoltaic devices based on shift currents. By analyzing 
effective models, here we outline simple design principles for the optimization of shift currents for 
frequencies near the band gap. Our method allows us to express the band edge shift current in 
terms of a few model parameters and to show it depends explicitly on wavefunctions in addition to 
standard band structure. We use our approach to identify two classes of shift current photovoltaics, 
ferroelectric polymer films and single-layer orthorhombic monochalcogenides such as GeS, which 
display the largest band edge responsivities reported so far. Moreover, exploring the parameter 
space of the tight binding models that describe them we find photoresponsivities that can exceed 
100 mA W“^. Our results illustrate the great potential of shift current photovoltaics to compete 
with conventional solar cells. 


Introduction - Cost-effective, high-performing solar 
cell technology is an essential piece of a sustainable en¬ 
ergy strategy. Exploring approaches to photo-current 
generation beyond conventional solar cells based on pn 
junctions is worthwhile given that their performance is 
in practice constrained by the Shockley-Queisser limit^. 
One of the most promising alternative sources of pho¬ 
tocurrent is the bulk photovoltaic effect (BPVE) or ‘shift 
current’ effect, a non-linear optical response that yields 
net photocurrent in materials with net polarization^“^°: 
Contrary to conventional pn junctions, the BPVE is able 
to generate an above band-gap photovoltage^^, poten¬ 
tially allowing the performance of BPVE-based photo¬ 
voltaics to surpass conventional ones. However, closed- 
circuit currents generated via the BPVE reported in the 
literature have typically been small compared to those 
generated in pn junction photovoltaics^^^^^. Recent in¬ 
terest in the BPVE also stems from the proposal that it 
may be at work in a promising class of materials for pho¬ 
tovoltaics known as hybrid perovskites^^, an extremely 
active field of research^®^^®. 

The fundamental requirement for a material to pro¬ 
duce a current via the BPVE is that it breaks inver¬ 
sion symmetry, allowing an asymmetric photoexcitation 
of carriers. But despite considerable case-by-case study 
of the BPVE, the necessary ingredients to optimize a 
BPVE-based solar cell are not sufficiently well under¬ 
stood. As with conventional solar cells, band gaps in 
the visible (1.1-3.1 eV) and large electronic densi¬ 
ties of states^^’^° are always beneficial. In addition, to 
produce a solar cell that responds to unpolarized sun¬ 
light, a highly anisotropic material must be used, since 
otherwise there is no preferred direction for the current 
to flow. But beyond these natural requirements, our only 
guiding knowledge is that the shift current depends ex¬ 
plicitly on the nature of the electronic wavefunctions^°’^^ 
and that it is not correlated with the material polariza¬ 
tion in any obvious way^'^ despite the fact that both shift 
currents and polarization originate from inversion sym¬ 
metry breaking. 



FIG. 1. Schematics of proposed shift current photovoltaics: a) 
3D structure of a solar cell built by stacking one-dimensional 
ferroelectric polymers, b) Simplified two-band tight binding 
model of a polymer, c) 3D structure of a solar cell made 
by stacking two-dimensional monolayers of a monochalco- 
genide. The inert spacers between layers prevent the restora¬ 
tion of bulk inversion symmetry, d) Simplified two-band tight- 
binding model for a monochalcogenide layer. 


In the current situation, a more generic understand¬ 
ing of what makes the BPVE strong is highly desirable. 
When tackling complex material science problems, strip¬ 
ping off all complications and optimizing the simplest 
model that captures the relevant physics often proves the 
best strategy, as shown for example in thermoelectricity 
studies^^^^^. In this work, we present simple design prin¬ 
ciples for BPVE optimization based on the study of an 
effective model for the band edges. With this model, 
the band edge shift current is given by the product of 
the joint density of states (JDOS) and a matrix element, 
both given by simple expressions in terms of a few model 
parameters. The simplicity of the model allows us to de¬ 
rive the main principle that band edges with semi-Dirac 
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type of Hamiltonians are the best starting point to ob¬ 
tain large band edge prefactors. In addition, by relating 
the effective model parameters to realistic tight-binding 
models, we can predict that several materials with the 
required band structure have larger shift currents than 
any reported so far. 

Results - In our search for materials we should look 
for large JDOS in systems where the band edge is closely 
aligned with the peak of the solar spectrum, around 1.5 
eV. Since the band edge always induces a Van Hove sin¬ 
gularity in the density of states, the requirement of a 
large peak in the photoresponse can be naturally better 
satisfied by low-dimensional materials, which generically 
present stronger singularities^^. Materials of one and two 
dimensions are therefore the focus of this work. Among 
one-dimensional materials, ferroelectric polymers are 
suitable candidates for shift-current photovoltaics: they 
strongly break inversion symmetry, some have suitable 
band gaps for photovoltaics applications^®^^®, and they 
can be produced in macroscopically oriented samples. 
For these reasons, we consider solar cells consisting of 
such polymer films, shown in Fig. 1(a). Two-dimensional 
materials^° also have great potential for photovoltaics, as 
shown by demonstration of a pn-junction photovoltaic ef¬ 
fect in dichalcogenide heterostructures'*^^^^, and in few- 
layer black phosphorus"^^. However, these well known 
2D semiconductors have vanishing shift currents because 
of either inversion or rotation symmetry. Group IV 
monochalcogenides have emerged in the past years as 
a new familiy of inversion-breaking, anisotropic 2D ma¬ 
terials with fascinating properties^®^^®, and interest is 
growing as thin films of all four members of the family, 
GeS50-53^ GeSe52.53^ SnS^hss and SnSe^^-ss, have now 
been isolated experimentally. In this work, we show that 
GeS is ideally suited to realize high values of the BPVE. 
Their GeS structure is shown in Fig. 1(c). 

To understand how to optimize the photoresponse, we 
first discuss how the shift current can be computed for 
a tight binding model, and then we proceed to apply 
this formalism to describe a generic band edge and the 
response of particular materials. 

Shift current - In this work we consider the shift 
current contribution to the BPVE and we shall use both 
terms interchangeably (note the BPVE can have other 
contributions as well®). With electric field Etj{uj) at fre¬ 
quency w and linearly-polarized in the b direction, the 
shift current is a DC response of the form® 

Ja = a^^\oj)Eh{cj)Eh{-oj). ( 1 ) 

Defining an intensity for each polarization, /o,b = 
ceo\Eb\‘^/2, we define the photoresponsivity as the 
current density generated per incident intensity Ja = 
which gives = 2a°'^^/ceo. Note that in con¬ 
ventional solar cells the current is also linear with inten¬ 
sity. For a D-dimensional system, takes the form^’® 

/ 7^ E - CO), ( 2 ) 


where C = Aggire^/fi^eoc, with c the speed of light, eo 
the vacuum permittivity, and gs = 2 accounts for the 
spin degeneracy. In what follows we set H = 1. Sum¬ 
mation of indices is explicitly indicated using the sum¬ 
mation symbol. The sum is over all Bloch bands, with 
= En — Em the energy difference between bands n 
and m and i^m = fn ~ fm the difference of Fermi occupa¬ 
tions, which we take at zero temperature. The integrand 
is 


-fnm = (3) 

where are the inter-band matrix elements of the posi¬ 
tion operator (or inter-band Berry connections), defined 
as r^m = * for n 7 ^ m and zero otherwise, where 

|n) is the eigenstate of band n. A semicolon denotes a 
generalized derivative r^m-a = 

where = i{n\dk^n) is the diagonal Berry connection 
for band n. 

Generic two band model - With the aim of de¬ 
scribing the shift current response of the band edge of 
a semiconductor, next we consider the shift current of a 
generic two band model. The Fourier transform of the 
real space Hamiltonian is performed with the choice of 
phases ^pk{x) = ^ En,! - R - x,) 

where 4>(x) is a localized orbital and is the position 
of site i in the unit cell. This choice is made in order to 
naturally incorporate the action of the position operator, 
see Refs.®®^®^. The Hamiltonian matrix takes the form 

+ E 


where cto is the identity matrix, ai = (Jx,cry,(7z are the 
Pauli matrices and eo and fi = fx, fy, fz are generic func¬ 
tions of momenta k (the momentum label is omitted to 
simplify notation). The conduction and valence bands 
are given by = eo -|- e, E 2 = cq — e, respectively and 
e = (Ei Note that this basis choice implies that 

the Hamiltonian matrix elements are not periodic in the 
Brillouin Zone, Hij (k -f G) 7 ^ Hij (k) with G a reciprocal 
lattice vector. 

To compute the shift current, the direct use of Eq. 3 
requires the evaluation of derivatives of Bloch functions, 
which can be difficult to compute numerically. Previ¬ 
ous works"^’^’® have addressed this problem with the use 
of identities that replace wavefunction derivatives with 
sums over all states of matrix elements of Hamiltonian 
derivatives. These identities are known as sum rules and 
rely on the fact that momentum and velocity operators 
are proportional in the plane wave basis p = mv, which 
is not true in the tight binding formalism. In this work 
we derived a generalized sum rule appropriate for tight 
binding models (see Methods section), from which the in¬ 
tegrand Eq. 3 can be evaluated for any two-band model 
in terms of the Hamiltonian derivatives only. The result 
is 


^12“ = - E - fmfi,bfj,a^-j)eijm, 

ijm 


( 5 ) 



3 




1.0 1.5 2.0 2.5 3.0 

cj(eW) 


FIG. 2. Frequency dependence of the components of photoresponsivity for different tight-binding models, computed from 
Eqs. 2 and 5: (a) Responsivity for a stack of disubstituted polyacetylene polymers with tight binding parameters ti = 2.85, 
t 2 = 2.15, A = 1.0 in eV, showing the square root divergence of the current at the band edges, (b) Various non-zero components 
of the responsivity tensor for a stack of 2D monochacogenides with parameters ti = —2.33, t 2 = 0.61, ts = 0.13, A = 0.41 in 
eV, and xo = 0.52A. A large peak is observed in at the band edge, (c) Responsivity for A = 0.8 eV, xo = O.OA, ts = 0 
and different hopping ratios |ti|/t 2 approaching the semi-Dirac limit. The emergence of a singularity is observed. In the three 
figures, solid lines show the shift current components as computed from the tight-binding model, and a dashed line in each 
subfigure shows the xxx shift current component as predicted by the effective low energy model valid near the edge, Eq. 9 


where the compact derivative notation fi^a = dk^fi and 
is used. Eq. 5 is one of the main results of this 
work. Several general principles to maximize the band 
edge shift current can be derived from this expression. 
A straightforward one is that, since this expression does 
not depend on eg, particle-hole asymetry does not influ¬ 
ence the shift current at all. Therefore eg is set to zero 
from now on. The additional term that appears only 
for tight binding models in this more general sum rule is 
fmfi,bfj,ab, wliich is absent in previous formulations. For 
a direct band gap, this term dominates the response ex¬ 
actly at the band edge, since to lowest order in k the first 
term always has constant contribution, while the second 
one is at least linear in k for any model due to the en¬ 
ergy derivative e^b- For this term to be finite, the three 
Pauli matrices in the Hamiltonian must have constant, 
linear, and quadratic coefficients, in any order. Satisfy¬ 
ing this low-energy constraint can be taken as another 
general principle in the search for materials with large 
shift current. 

More explicit guidelines can be obtained by consider¬ 
ing an explicit low-energy model with a direct band gap 
at a time reversal invariant momentum. Expanding the 
Hamiltonian around it we get 

H -“t“ G-xk^ “t“ Gyky -\- Gxykxky'^CTx 

-\-Vpkx(Xy (A (^x^x Pxy^x^y)^z- (6) 

Time reversal symmetry iJ*(—k) = H(k) prevents 
quadratic terms in ay, and we have taken the linear term 
to be in the x direction without loss of generality. Note 
this type of linear term requires the breaking of any Cn 
rotation symmetry with n > 2. The band gap of this 
model is Eg = 2ek=g. Evaluating 5 we get 


= -^{axA - px5) + 0(y), (7) 

AT(^) = - Pxyd) + 0{k^), (8) 

while 1^2^ = = 0 + O{k^). Also note that in order to 

have a non-zero shift current quadratic terms in ax or a^ 
are required. In 2D, the fact that is in general non¬ 
zero means that the current need not be in the direction 
of the electric field polarization. 

The shift current close to the band edge can now be 
obtained by substituting Eqs. 7-8 into Eq. 2, which gives 

E'’\u:) = C (tc-Eg)/Eg<l (9) 

where N{uj) = J dk^ S{uji 2 — uj)/{2Tr)^ is the JDOS. 
Eq. 9 provides an analytical formula for w close to the 
band edge for a very general class of models. This simple 
expression allows one to disentangle the contributions of 
the shift current integrand and the JDOS and hence to 
optimize them independently. 

To maximize the response we therefore require band 
structures where the JDOS has a strong singularity. It 
is well known that in the ID case, the generic JDOS 
diverges as a square root, N{u}) cx (w — Eg)~^/‘^. ID 
systems such as polymers or nanowires or systems in the 
quasi ID limit will in general have a large response. In 
2D, the band edge JDOS has a finite jump of N{uj) = 
[mx’oiy')^^'^I2 tx , where rai are the average effective masses 
for valence and conduction bands. A singular N{lo) thus 
occurs in 2D when the inverse effective mass vanishes. In 
the effective model in Eq. 6, this happens when <5 = 0, 
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which realizes what we may call a gapped semi-Dirac 
dispersion®^, since the coefficients of Uy and Ux are linear 
and quadratic in momentum, respectively. In such a case 
we have N{ijj) oc (w — (full expressions for N{uj) 

may be found in the Methods section). 

For materials with large JDOS, the current can be fur¬ 
ther enhanced by appropriately tuning the parameters in 
Eqs. 7-8. This is most easily discussed if these parame¬ 
ters can be related to microscopic lattice models. In the 
next section, we discuss tight-binding models for simple 
materials that realize the described types of band struc¬ 
tures. 

Material realizations and lattice models - As a 

realization of the ID case, we consider ferroelectric poly¬ 
mers that break inversion symmetry such as polyvinyli- 
dene fluoride or disubstituted polyacetilene®®’®®’®®. This 
system is described by the tight-binding model schemat¬ 
ically shown in Fig. 1(b), defined in terms of two types 
of hoppings, ti and t 2 , alternating on-site potentials ±A, 
and orbital centers at a; = 0 and x = xq. With our 
choice of basis functions, the Hamiltonian is specified by 
fx + ify = —and fz = A, where 
a = 10 A is the lattice constant and the distance be¬ 
tween closest neighbors is®® a;o = 0.48a. For estimates of 
the tight binding parameters, we consider the example of 
disubstituted polyacetilene that was experimentally re¬ 
alized in Ref. 38, with a band gap of 2.5 eV. For regu¬ 
lar polyacetilene, where A = 0, the hopping parameters 
and band gap have been estimated as®® H = 2.85 eV, 
t 2 = 2.15 eV, Eg = 1.4 eV. Assuming the same hop¬ 
ping for the disubstituted version, we use A = 1.0 eV to 
match the observed band gap. Note that the dispersion 
does not depend on xq. 

Using Eq. 2 and Eq. 5 we can now compute the shift 
current for this ID model. Expanding about the low en¬ 
ergy momentum kx = 7r/a and performing a constant 
rotation of the Pauli matrices, we obtain an effective 
model as Eq. 6 with parameters fcy = 0 and 6 = ti — t 2 , 
VF = (U - ^ 2 ) 3^0 + t2a, ax = [t2ia - xq)^ - tixl]/2. 

To be able to compare the responsivity of these ma¬ 
terials to that of a 3D system, we consider a stack of 
polymers as depicted in Fig. 1(a), separated by a dis¬ 
tance d which we take to be equal to the lattice constant 
of the polymer d = a. The photoresponsivity is then 
The typical photoresponsivity spectrum 
of this model with this convention is shown in Fig. 2(a). 

For the 2D case, we require a layered material that 
breaks both inversion and rotational symmetries. The 
most popular of the recently isolated 2D semiconduc¬ 
tors break either inversion (BN, M 0 S 2 ) or rotational 
symmetries (black phosphorus®^, ReS' 2 ®®), but not both. 
An inversion symmetry breaking version of the strongly 
anisotropic black phosphorus, a group V element, can be 
obtained combining elements of the IV and VI groups. 
These group IV monochalcogenides, such as GeS, are 
predicted to be stable in the monolayer form with the 
orthorhombic structure of black phosphorus"*®’^®. 

These materials can be described with a tight binding 


model similar to the one used for black phosphorus®®^®®. 
While the GeS unit cell contains two Ge-S pairs at dif¬ 
ferent heights, a unit cell with a single Ge-S pair can 
be used when the physics to be probed is insensitive to 
the heights of the atoms (see Methods for a detailed 
explanation). The two band Hamiltonian is specified 
by fx + ify = -b t 2 ^'(k) + < 3 $*(k)], where 

Xq = (xojO) and ‘I’(k) = -band fz = A. ai 

and a 2 are the lattice vectors. See Fig. 1(d) for the defini¬ 
tion of the hopping integrals. Again note the dispersion 
is independent of Xg- The specific values of the tight- 
binding parameters for GeS have been obtained by fit¬ 
ting an ab-initio calculation as described in the Methods 
section, where the coefficients of the low energy model 
near the band edge are also shown. Note in this lattice 
structure there is a mirror symmetry y —)■ —y, which is 
represented as the identity, and restricts axy = Pxy = 0. 
(This is so because both conduction and valence bands 
are even under the symmetry, as it also happens in black 
phosphorus. This is also the result of our ab-initio calcu¬ 
lation.) This symmetry still allows a linear term of the 
form kxCTy, crucial for the semi-Dirac type of band struc¬ 
ture. In this model, the semi-Dirac limit is realized when 
ti = — 2(^2 + 

We consider a stack of monolayers separated by d = 
a, as shown in Fig. 1(c). In this case, we consider an 
inert spacer layer between the GeS layers to avoid the 
restoration of inversion symmetry that would occur if we 
were to stack GeS into its natural bulk form. The 3D 
photoresponsivity of this model, given by = K 2 I) /d, 
is computed using Eqs. 2 and 5. To make contact with 
the ID case we consider a stacking distance d = a = 
(|ai| -b |a 2 |)*/^ and xg = 0.18a. The results are shown in 
Fig. 2(b). We see that both and are in general 
finite, and the polarization average is also finite due to 
the strong anisotropy. 

The response of the monochalcogenides is large be¬ 
cause they are close in parameter space to the gapped 
semi-Dirac Hamiltonian. This is best illustrated by con¬ 
sidering the evolution of a fictitious system where the 
hoppings are tuned (with tg = 0 for simplicity) to the 
semi-Dirac case |ti|/t 2 = 2, where the divergence of the 
response is clearly appreciated. This evolution is shown 
in Fig. 2(c). 

Further optimization - After describing the repre¬ 
sentative tight-binding models with large JDOS, we may 
now address a more systematic analysis of the photore¬ 
sponsivity. First, we consider exploring the phase di¬ 
agram of the monochalcogenides by sweeping |ti|,t 2 in 
parameter space while the band gap is fixed at 1.89 
eV by choosing A appropriately and ts = 0 for sim¬ 
plicity. Fig. 3(a) shows the polarization averaged pho¬ 
toresponsivity, Kx = (k“®^ -b K^'^y)/2, for the parameters 
xg = 0.18a and 0 = 0.69. This phase diagram summa¬ 
rizes nicely the most physically relevant regimes where 
the shift current is large due to a divergent JDOS, namely 
the ID dimensional limit where |ti| ^ t 2 , and the semi- 
Dirac regime where |ti| ~ 2t2- In this phase diagram. 
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FIG. 3. Phase diagrams for monochalcogenide layer tight binding model: (a) Polarization-averaged photoresponsivity in the 
^-direction, at the band gap frequency plotted as a function of hopping parameters |ti| and t 2 , keeping the band gap fixed 
at 1.89 eV by tuning A accordingly. The Ge-S distance is xo = 0.52 A and ta = 0. The location of GeS on the phase diagram 
is marked by a white circle with blue outline. Regions for which the gap cannot be kept at 1.89 eV are left white, (i) and 
(ii) show bond strengths in the limits where |ti| 3> t 2 and |ti| <C t 2 , respectively, to illustrate the two extremes of the phase 
diagram, (b) Polarization-averaged photoresponsivity in the i-direction, R^, at the band gap frequency plotted as a function 
of the Ge-S distance xq in units of a = (af -|- and ratio of hopping parameters |ti|/t 2 . Here, A and t 2 are set to GeS 

values of 1.1 eV and 0.61 eV, respectively. The location of GeS on the phase diagram is marked by a white circle with blue 
outline, (iii) and (iv) show two extreme cases of the phase diagram, where xo is large and small, respectively. 


the point corresponding to ti and t 2 of GeS is shown as 
a white circle with blue outline. 

Next we illustrate a very important feature of the be¬ 
havior of the shift current integrand. Eqs. 7-8 depend 
generically on the hoppings and lattice parameters. The 
energy does not depend on the parameter xq, but the 
wavefunctions do. In Fig. 3(b), we show the peak pho¬ 
toresponsivity as a function of |ti|/t 2 and xq. A large 
response is observed in the semi-Dirac limit |ti|/t 2 2. 

However, a very strong dependence on Xq and even a 
sign change is also observed. The dependence on xq dra¬ 
matically illustrates the fact that the shift current de¬ 
pends not only on the band structure, but also on the 
wavefunctions. This can be seen explicitly in the fact 
that the effective mass m~^ = 4a^tit2/7?g is indepen¬ 
dent of Xq, but the combination vpa^ appearing in the 
shift current integrand is not. In particular ax vanishes 
for Xq = ax/Y + (|G/ 2 t 2 |)^^^], which means that regard¬ 
less of the JDOS, the band edge response can actually 
be zero. This behavior is characteristic of Berry con¬ 
nections, which depend explicitly on the positions of the 
sites in the unit cell. 

Discussion - In this work, we have shown how an 
effective model for the band edge enables a clean sepa¬ 
ration of the two factors that contribute to a large shift 
current: the standard JDOS and the shift current matrix 


element. This model also allows us to readily identify ma¬ 
terials with semi-Dirac-like Hamiltonians as those where 
both factors can be made large. Several other general 
conclusions can be drawn from the form of the effective 
shift current integrand in Eqs. 7-8. First, since the l/w^ 
factor becomes l/£^g at the band edge, materials with 
smaller gaps are expected to have larger shift currents. 
A second conclusion is that while looking for materials 
with large JDOS is a good guiding principle, the shift cur¬ 
rent integrand depends on other microscopic details that 
can change the response dramatically. Within our simple 
model, the shift current can be maximized by bringing 
the two sites of the unit cell closer together, which is 
a requirement that the monochalcogenides satisfy well. 
Materials that may perform even better than GeS may 
be searched for exploring different chemical compositions, 
alloying, or by strain engineering. 

Our results were made possible by the derivation of 
a new sum rule appropriate for tight-binding models. 
With this sum rule, our work can be easily extended 
to tight-binding models with more than two bands, or 
systems where the minimum direct gap is not at a time- 
reversal invariant momentum. We expect that the for¬ 
malism developed here will provide the necessary link to 
combine ab-initio methods with effective models, allow¬ 
ing for more in-depth, systematic study of shift current 
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photovoltaics. 

Our results should be compared to known ferroelectric 
materials that have been recently studied. In the visible 
range of frequencies, w < 3 eV, we find peak values of 0.1 
mAW“^ in BiFeOs 1 mAW“^ in hybrid perovskites^^ 
and a maximum 10 mAW“ ^ in BaTiOs^^ or NaAsSe 2 ^^. 
The realistic materials that we propose present larger 
responsivities, with the additional advantage that the 
peak is by construction at the band edge. Moreover, 
as Fig. Fig. 2(c) and Fig. 3(b) show, peak responses 
on the order of several hundreds of mAW“^ could be 
achieved with materials closer to the semi-Dirac regime. 
To compare with conventional photovoltaic mechanisms, 
the total current per intensity of a crystalline Si solar cell 
exposed to sunlight is about 400 mAW“^^°. 

Given these numbers, our work is a sign that shift cur¬ 
rent photovoltaics capable of surpassing conventional so¬ 
lar cells may be close at hand, and a push to investigate 
their full potential using methods discussed in this work 
- along with established techniques - is warranted. We 
believe that the simple principles derived in our work will 
serve as a guide for both theory and experiment in the 
development and optimization of the next generation of 
shift current photovoltaics. 

Methods 

Shift current - To make contact with previous work, 
we note the shift current integrand in Eq. 3 is sometimes 
expressed in terms of the phase of the inter-band matrix 
element as = Ir^rn^R^t where 


= + ( 10 ) 

is known as the shift vector. The response to a natu¬ 
ral light source such as sunlight, which is unpolarized, 
is obtained by averaging over polarization. Taking 
E{0) = |i?|(cos0,sin0) we have 

Ja= J ^Ja = = nJo- ( 11 ) 

Sum rule - The expression for the shift current pre¬ 
sented in the main text can be obtained by the use of a 
sum rule for the quantity which is obtained from 

the identity 


dkbdka {n\H\m) = Snmdkbdk^En- (12) 

Evaluating both sides explicitly for n ^ m, the identity 
can be expressed as 


^nm;b 


nm nm 


yb 

nm nm 


— W. 


ab 


E( 

p^n,m 


UJn 


yb ya 

^np ^pm > 
^nv 


n ^ m (13) 


where = {n\dkbH\m) are the velocity matrix ele¬ 
ments, = {n\dkbdk^H\m) and 

ujnm = En — E^- lu the evaluation, we used 


(t^ 

\' mnJ nmf 

'^nn ~ dk^En, 


(14) 

(15) 

n ^ m (16) 


The first equality follows from dk (n|m) = 0 if m n, 
while the last two follow from dk^ {n\H\m) = SnmdkaEn- 
Note this sum rule contains the extra term com¬ 

pared to Ref.®, where E[ = p®/2m -I- V{x) and = 
Im which has no off diagonal component. Quite 
importantly, the term in tight binding models is the 
one responsible for all band edge contributions. Also note 
that it has been argued before that = 0 for a two 
band model^, which is actually only true if = 0. 

Two band model - For the case of two bands, m = 
1, n = 2 the use of the sum rule for the shift current 
integrand in Eq. 3 leads to the simplified expression 


jabb 

^nm 


-Im 


UJ 


12 


-V^tV 


21 ^ 12(^11 
2e 


b ] 
22 ) 


+ Vn^W 


ba 


(17) 


To evaluate this expression we compute the wave func¬ 
tions of E[ 

V’n = ;^(-r7\/e-»7/^ ,e*'^''\/e + »7/z), (18) 


with n = 1,2, ry = (—1)”, and (()k = &rctan{fy/fx). The 
required matrix elements are 


^^21 = (V'2|(eo.a2:+ alV'l) 

i i 

W 12 = {iplliio,ab^ + ^crifi,ab\'ip2) = ^ ( 20 ) 


where the off diagonal matrix element Si = {tpi\o'i\'ip 2 ) is 


Si = (^ cos (fiu + i sin ^-f sin ((ik 


— l cos 0k) ~ 




( 21 ) 


and the diagonal velocity matrix elements are computed 
from Eq. 15. The imaginary part in Eq. 17 can be taken 
using Im [s*Sj] = — 'Ylim ^ijmfm/^ and this leads to Eq. 5 
in the main text. 

Joint density of states - To compute the JDOS, we 
first start with the ID case. Close to the band edge, we 
expand the energies of conduction and valence bands as 
Ei « Ei{0) + k‘^l2mi^x, so that u;i 2 = Ei — E 2 ~ E^ + 
k\l2mx where the total effective mass m~^ = -I- 

\m 2 ,x\~^ is given by 

= 4(u| + 2a,,5 + 2j3x^)/E^, (22) 


and solve for k{(jj) = yj2mx{oJ — Eg). Rescaling 2mx we 
get 
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N 


ID 


, . f- - f dk 6(k I 

(a;) = V2^y - — 


dk 5{k ± k{uj)) 

y/2mx d{uj — Eg) 

27T 


(23) 


where we get the expected ID singularity. For the generic 
2D case, again we expand wi 2 ^ Eg + k^/2mx + ky/2my, 
where rux is still given by Eq. 22 and 

m-^ =8{ayS + pyA)/Eg, (24) 

We consider the case when > 0, my > 0, so that the 
minimum does lie at k = 0. By rescaling 2mx and 2my 
we get in polar coordinates 


N 2 D 



kdkdO S{k — k{oj)) 

(27r)^ \2k\ 


y/m^my 


9[ijj 


Es), 


(25) 


which is the expected constant result. Finally, the semi- 
Dirac case occurs in 2D when = 0, which in the 
absence of second neighbor hopping occurs exactly at d = 
0. In this case, we keep the complete expression for a;i 2 = 
{{a^kl -I- ctyky)'^ + Vpk'l -|- A^)^/^. In polar coordinates 
we have 


constructed in a standard way by projecting into hydro- 
genic s-like and p-like orbitals on both Ge and S atoms 
along with two s-like orbitals in the vacuum region that 
are needed to represent the vacuum states. The frozen 
window for the disentanglement procedure spans up to 
6.2 eV above the Fermi level. The crystal structure of 
GeS is orthorombic with space group Pnma (No. 62) 
and lattice vectors li = (?i,0) and I 2 = ( 0 ,/ 2 ), with 
li = 4.53 A and I 2 = 3.63 A and contains two Ge and two 
S atoms. The structure can be seen as two GeS zigzag 
chains separated by a height oih = 2.32 A. The ab-initio 
results for the conduction and valence bands near the F 
point are shown in Fig. 4 and have mostly Pz character. 



2n In 


atSd f kdkdO 5{k — k{w)) 

^ J (2^)2 |afcCCi2| 

We now rescale Oa,, ay instead, solve for k 


FIG. 4. Tight binding fit to ab initio for GeS: Dispersion of 
conduction and valence bands of GeS near F computed ab- 
initio (red dots). A black line shows the tight binding ht for 
comparison. 


k{uj) = [-vj./axCos^0±{vj,/alcos'^9 + uj^-El)E‘^]/2, 
and get 

j^SD ^ f dkd9 S{k — fc(cLi)) 

A^ya^ay J ( 211 )"^ {vp/a^ cos'^ 9+ uj'^ — Ej)'^A 

^ _£(i)_ a;6l(a; - Eg) 

4r(|)(27r)3/2|aa,|ya^UF (w^ - E2)i/4 

Ab-initio calculation and tight binding fit for 

GeS - Due to the lack of tight binding models for 
monochalcogenide materials'^^^e^ have derived the 
tight binding parameters by fitting the electronic struc¬ 
ture of GeS ab-initio. We used the PBE^^ approxi¬ 
mation to the exchange correlation functional, ultrasoft 
pseudopotentials ,^2 Quantum-ESPRESSO^^ and Wan- 
nier90^^ computer packages. The cutoff for electron 
wavefunction is set to 40 Ry and cutoff for electron den¬ 
sity to 200 Ry. Internal coordinates and in-plane lattice 
constants were fully relaxed. Vacuum region between 
repeating images of GeS monolayers is 17 A. Wannier 
functions were constructed from a 12x12 regular k-mesh 
grid. The maximally localized Wannier functions were 


This system can be effectively described with a two site 
tight binding model. This can be done because the lattice 
structure has glide symmetries with mirror reflection 2 ; —>■ 
—2 and translations Si = {ax, ay) and 02 = {ax,—ay), 
with Qx = l\/2 and ay = ^ 2 / 2 . When the out of plane 
positions of the atoms are not relevant for the problem of 
interest, one can define a smaller two site unit cell where 
the glides play the role of lattice vectors (as it is done in 
black phosphorus®®). The Ge and S sites in this effective 
tight binding model are located at (0, 0) and {xg, 0), with 
Xg = 0.62 A. This is the tight binding model employed in 
the main text. The parameters of this model are obtained 
from the ab-initio calculation as follows. 

Since our aim is to model faithfully only the low energy 
bands around the Gamma point, it will suffice to consider 
a single Pz orbital per site in the tight binding model. 
The minimal model parameters are the on-site potential 
difference A between Ge and S Pz orbitals and the three 
nearest neighbors hoppings ti, with f = 1, 2, 3, which are 
all between Ge and S atoms. In addition, to reproduce 
the small particle-hole asymmetry of the gap, we also 
consider two further neighbor hoppings t'l and <2 which 
connect Ge-Ge or S-S pairs (we assume the same values 
for both species to simplify). 


















The tight binding Hamiltonian takes the form H = 
eo + SitTi/i(k) with coefficients 


The Berry curvature around T for the tight binding 
model is given by 


eo = — 2t'i (cos ai • k + cos el 2 ■ k) 

— 2^2 cos(ai — a 2 ) • k, (28) 

/. + ifv=- + t 2 $(k) + t3$*(k)], (29) 

A =A, (30) 

where, as defined in the text, $(k) = Our 

tight binding fit is intended to reproduce faithfully the 
bands and wavefunctions close to the band edge, where 
the effective low energy model applies. This model is 
given by 


_ VF{ay^ - f 3 y 6 ) 

“ (A^ ' ^2^3/2 


(40) 


. ^ 2 ) 3/2 

Since H vanishes at the origin, we take dyO, as one extra 
input for the ht. The quantum metric is defined as 


AyAjy, (41) 

The only non-vanishing component of the quantum met¬ 
ric at fc = 0 is given by 


Qxx — 


4(A2-hd2)’ 


(42) 


H = ilxkl -I- lykl)I -I- (d -I- aa;kl + aykl)a^ 

+ vpkxcry + (31) 

where a constant term is omitted as it can be absorbed in 
the chemical potential. The effective model parameters 
are related to the tight binding parameters as 


lx — (32) 

^y = {2t\+AQal, (33) 

S = ti — 2t2 — 2t3, (34) 

vf = -2a^{t2 - h) - (ti - 2t2 - 2t3)xo, (35) 

otx = t2{ax - xof - tixl/2 + tsiax + xq)^, (36) 

ay = {t2+t3)al. (37) 


The key to obtain a reliable tight binding parametriza- 
tion is that, since the shift current depends sensitively on 
the actual wavefunctions, the tight binding model should 
be fitted to wavefunction dependent quantities in addi¬ 
tion to the band energies. The simplest gauge invari¬ 
ant quantity that depends on wavefunction phases is the 
bracket of two covariant derivatives 

Qflu — {k^fi^k\^u'^k} 5 (^^) 

with Dy, = dy, - iAy, with Ay = i {uk\dyUk) the Berry 
connection. The real and imaginary parts of this tensor 
are known as the Berry curvature and the quantum met¬ 
ric. A Ht that reproduces this tensor correctly in addition 
to band energies ensures that the wavefunction structure 
around the T point is correctly accounted for, so that any 
other gauge invariant quantity computed in the effective 
model should be the same as that computed ab-initio. 


The Berry curvature H(fc) is defined as 

n{k) = eyAkn[{9fiUk\duUk)] = V x A. (39) 


so we take as another extra input for the fit. 

In summary, we take as ab-initio input parameters the 
gap, the four effective masses, and the lowest order Berry 
curvature and quantum metric, dyD, and gxx- The dif¬ 
ference in effective masses for electron and hole bands, 
accounted for the term eo, can be fitted independently 
with the hoppings t'l and Since eg has no impact in 
the shift current response, the hoppings t'l and ^re not 
considered in the main text. The rest of the input is fit¬ 
ted with ti, t 2 and ts, the on-site potential A and xq, 
and the results of the fit are shown in Table I. While xg 
is in fact known from the lattice structure of GeS to be 
0.62A, obtaining it independently from the tight bind¬ 
ing fit, which gives a close value of 0.52A provides an 
additional check of the validity of the model. 
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1 Ab-initio input parameters | 



'^X,C 


'^T^y,c 

dyQ, 

Qxx 

1.89 eV 

-0.064 eV^A-2 

0.079 eViA-2 

-0.340 eV"iA"2 

0.171 eV"2A"2 

3.565 A® 

2.529 A2 

1 Tight binding parameters | 

A 


t2 

ts 

t'l 

t2 

Xo 

0.41 eV 

-2.33 eV 

0.61 eV 

0.13 eV 

0.07 eV 

-0.09 eV 

0.52 A 


TABLE I. Table of ab-initio and tight binding parameters for GeS: First row: input ab-initio parameters. Second row: Tight 
binding parameters obtained from the fitting. 
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